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ABSTRACT 

We perform a stability test of triaxial models in Modified Newtonian Dynamics 
(MOND) using N-body simulations. The triaxial models considered here have den- 
sities that vary with in the center and r~'* at large radii. The total mass of 
the model varies from IO^Mq to IO^'^Mq, representing the mass scale of dwarfs to 
medium-mass elliptical galaxies, respectively, from deep MOND to quasi-Newtonian 
gravity. We build triaxial galaxy models using the Schwarzschild technique, and evolve 
the systems for 200 Keplerian dynamical times (at the typical length scale of 1.0 kpc). 
We find that the systems are virial overheating, and in quasi-equilibrium with the re- 
laxation taking approximately 5 Keplerian dynamical times (1.0 kpc). For all systems, 
the change of the inertial (kinetic) energy is less than 10% (20%) after relaxation. 
However, the central profile of the model is flattened during the relaxation and the 
(overall) axis ratios change by roughly 10% within 200 Keplerian dynamical times (at 
l.Okpc) in our simulations. We further find that the systems are stable once they reach 
the equilibrium state. 

Key words: galaxies: kinematics and dynamics- methods: N-body simulations 
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1 INTRODUCTION 

Elliptical galaxies are often triaxial and appear stable. A 
triaxial equilibrium is non-trivial to build dynamically espe- 
cially for a system with a cuspy profile of the light and/or the 
dark halo. The main objective of this work is to test whether 
triaxial models of galaxies are stable in Modified Newtonian 
dynamics (MOND, Milgrom 1983). Extensive studies about 
the stability of triaxial models have been performed in stan- 
dard Newtonian gravity (see below), however, there is no 
literature on this topic in MOND. 

For Newtonian physics, it has been three decades of 
studies on constructing a self-consistent model for triaxial 
galaxies since Schwarzschild numerically presented the tri- 
axial Hubble profile in 1979 (Schwarzschild 1979; 1982). De- 
spite its original application to a modified Hubble profile, 
the method of Schwarzschild (1979) is still widely used for 
testing the self-consistency of various models for the density 
distribution in galaxies. For instance, Statler (1987) showed 
that the perfect triaxial Kuzmin (1973) profile and the de 
Zeeuw & Lynden-Bell (1985) profile are also self-consistent. 
Those models have constant density cores, however, observa- 
tions showed that elliptical galaxies have non-constant cores 
(MoUer, StiaveUi, & Zeilinger 1995; Crane et al. 1993; Jaffe et 



al. 1994; Ferrarese et al. 1994, Lauer et al. 1995), i.e. the sur- 
face brightness increases quickly towards the central region 
of the galaxies. Almost all elliptical galaxies have power-law 
cusps p ~ r ' with 7 ranging from 1 to 2 for High Surface 
Brightness to Low Surface Brightness elliptical galaxies in 
the central region. Spherical models with a fixed value of 
7 have been proposed, e.g. a 7 = 2 model by Jaffe (1983) 
and a 7 = 1 model by Hernquist (1990). Today such models 
are rather discussed within a family of density distributions 
with 7 being a free parameter (Dehnen 1993, Carollo 1993 
and Tremaine et al. 1994). In this regard, Merritt & Fridman 
(1996) tested the modified Dehnen profile. 



p{r) = 



(3 - 7)M 



47ra6c rT(l + r)' 



— ,0<7<3, 



(1) 



2 + (f )2, (c < fe < a), a, b and c 
and short axis of the ellipsoids. 



where r = ^(f )2 + (|, 
is the long, intermediate 
They found that triaxial galaxies with central density cusps 
(7=1) were in equilibrium and self-consistent in Newto- 
nian dynamics. The subsequent work by Capuzzo-Dolcetta 
et al. (2007) proved that a two-component triaxial Hernquist 
system, including a baryonic component plus a Cold Dark 
Matter (CDM) halo are also self-consistent. 
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Modified Newtonian Dynamics (MOND) - proposed by 
Milgrom (1983a,b) as an alternative gravity theory - was 
initially designed to abandon the need for that yet-to-be- 
discovered dark matter that (possibly) accounts for as much 
as 85% of all matter in the Universe. MOND, on the other 
hand, perfectly predicts the rotation curves of galaxies as 
well as the TuUy-Fisher relation in the absence of CDM 
(McGaugh et al. 2000; McGaugh, 2005). Indeed, MOND 
successfully matches the observations on a wide range of 
scales, from globular clusters (Angus & McGaugh 2008, in 
preparation) to different types of galaxies including dwarfs 
and giants, spirals and ellipticals (Milgrom 2007; Gentile 
et al. 2007; Milgrom & Sanders 2007; Famaey & Binney 
2005; Sanders & Noordermeer 2007; Angus 2008a). The de- 
velopment of several frameworks for a relativistic formula- 
tion of MOND (Bekenstein 2004; Sanders 2005; Bruneton & 
Esposito-Farese 2007; Zhao 2007; Skordis 2008) enabled the 
study of the Cosmic Microwave Background (CMB) (Skordis 
et al. 2006; Li et al. 2008), cosmological structure formation 
(Halle & Zhao 2008; Skordis 2008), strong gravitational lens- 
ing of galaxies (Zhao et al 2006; Chen & Zhao 2006, Shan 
et al. 2008) and weak lensing of clusters of galaxies (Angus 
2007; Famaey et al. 2007a). As a dynamically selected ref- 
erence frame, external fields break the Strong Equivalence 
Principle (Bekenstein & Milgrom 1984; Zhao & Tian 2006; 
Famaey, Bruneton & Zhao 2007b, Feix et al. 2008a,b). Con- 
sequently, the rotation curve, escape speed and morphology 
of galaxies are determined by both the background and the 
internal gravity (Famaey et al. 2007b; Wu et al. 2007, 2008). 
Despite its great success we need to accept that even MOND 
cannot do well without dark matter completely: a recent 
study utilizing a combination of strong and weak lensing 
by galaxy clusters showed that MOND requires neutrinos of 
mass 5 — 7eV (Natarajan & Zhao 2008). And to be consistent 
with (dark) matter estimates of galaxy clusters and obser- 
vataions of the CMB anisotropic spectrum (as well as the 
matter power spectrum), MOND requires neutrino masses 
of up to lleV (Angus 2008b). One theory capable of accom- 
modating both these requirements is that of a mass- varying 
neutrino by Zhao (2008). 

In this paper, we utilize a numerical solver for the MON- 
Dian analog to Poisson's equation to study the stability 
of triaxial galaxies in MOND. The code named NMODY 
has been widely applied to different problems: it has been 
applied to study dissipationless collapses, showing that 
the end-products are consistent with several observations 
(Nipoti, Londrillo & Ciotti 2006; Nipotti, Londrillo & Ciotti 
2007a). The code has also been used to study various impor- 
tant aspects of galaxy formation. Nipoti et al. (2007b) and 
Ciotti et al. (2007a, b) found that phase mixing is less effec- 
tive and the timescale of galaxy mergers is longer for MOND 
than for CDM. Recently, Jordi et al. (2009) and Haghi et al. 
(2009) applied the external fields into the NMODY code and 
studied the internal dynamics of distant star clusters. Fur- 
ther, MOND also produces stronger bars than CDM (Tiret 
& Combes 2007), and hydrodynamical simulations of spheri- 
cal bulges indicated that there are tight correlations between 
bulge mass, central black hole and stellar velocity dispersion 
in MOND (Zhao et al. 2008). These differences and simi- 
larities to CDM simulations immediately lead to the ques- 
tion of the stability of triaxial systems in MOND as realistic 
galaxies are not spherically symmetric objects. Wang et al. 



(2008) recently found that the self-consistency of a triaxial 
cuspy centre 7 = 1 also exists for MOND. By extending 
the original Schwarzschild method and weighting the orbits 
during the generation of the Initial Conditions (ICs) for N- 
body simulations it is possible to study the stability and 
future evolution of these density models (Zhao 1996). This 
method proved successful in, for instance, creating equilib- 
rium ICs for a fast-rotating, triaxial, double-exponential bar 
reminiscent of a steady-state Galactic bar (Zhao 1996) when 
evolved forward in time using a Self-Consistent-Field code 
(Hernquist & Ostriker 1992). 

Whether there are stable galaxy models in MOND is a 
lacuna in the studies. It is important to build stable galaxy 
models for dynamical studies. Here we will expand upon pre- 
vious work by studying the stability and evolution utilizing 
direct N-body simulations. Our target of study will be an iso- 
lated triaxial galaxy with a mild cusp of 7 = 1 in the centre 
within the Bekenstein- Milgrom MOND theory (1984). We 
use the same density models applied in Wang et al. (2008), 
with total mass ranging from 10^" Af© to 10* M0, respec- 
tively, representing medium-mass elliptical galaxies down 
to dwarf ellipsoidals, which are in quasi-Newtonian to deep 
MONDian gravity. We generate the ICs utilizing the method 
outlined in Zhao (1996) and our N-body simulations con- 
firm that these systems are (initially) in quasi-equilibrium 
and relax on a rather short time scale of only a few Keple- 
rian dynamical times (1.0 kpc) (see below, sub-section 13. ip . 
The systems quickly reach a state of equilibrium, consistent 
with the results of Wang et al. (2008). The inertial energy 
changes by less than 10% and the kinetic energy by less 
than 20% during the relaxation process. At the same time, 
the initial 7 = 1 cusps are flattened. After the relaxation, 
the systems remain stable. We further note that the triaxi- 
alities of the systems do not change significantly during 200 
Keplerian dynamical times (1.0 kpc). Moreover, the scalar 
Virial theorem is valid at any time. 



2 MODELS, SCHWARZSCHILD TECHNIQUE, 
AND ICS FOR N-BODY 

2.1 Poisson's equation in MONDian 

The MONDian Poisson's equation can be written as 
(Bekenstein- Milgrom 1984): 



|V$j 
fflo 



(2) 



where <I> is the MONDian potential generated by the mat- 
ter density p. For the gravity acceleration constant, we use 
ao = 3600 km^s~^kpc~^, which is same as adopted by Mil- 
grom (1983a,b), Sanders & McGaugh (2002) and Beken- 
stein (2006). The so-called MONDian interpolation func- 
tion jj, is approaching 1 for |V$| >> ao (Newtonian limit) 
and fi for |V<I>i « ao (deep MOND regime), and 

the gravity acceleration is then given by y^oogjv, taking the 
place of the Newtonian acceleration gjv = V^jv at the same 
limit. For our simulations we chose the 'simple' /^-function in 
the form of (Famaey & Binney 2005; Zhao & Famaey 2006; 
Sanders & Noordermeer 2007) 



li{x) = 



1 + a:: 



(3) 
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Furthermore, we use the density distribution given by Eq. [TJ 
choosing 7 = 1, and M being the total mass ol the system. 
For our simulations, we choose the ratios a : 6 : c = 1 : 0.86 : 
0.7, with a = Ikpc. 

2.2 Initial Potential 

A very important step in our calculation is to solve Pois- 
son's equation in MOND (cf. Eq. [2]). This is achieved via 
numerical integration utilizing the N-body code NMODY 
(Ciotti et al. 2006; Nipoti et al. 2007a) on a spherical grid 
of coordinates {r,0,tp). To this extent we applied a grid of 
256 X 64 X 128 cells. Note that we yet do not evolve the sys- 
tem forward in time; we simply extract the potential of our 
(static) density distribution and use it for the Schwarzschild 
method detailed below. 



2.3 Schwarzschild technique 

Since Schwarzschild (1979, 1982) pioneered the orbit- 
superposition method to construct self-consistent models of 
galaxies, this technique has been widely applied in dynami- 
cal studies (e.g., Zhao 1996; Rix et al. 1997; van der Marel et 
al. 1998; Kuijken 2004; Binney 2005; Capuzzo-Dolcetta et al. 
2007). The essence of this method is to sample phase-space 
with a large number of orbits. Properly assigning weights to 
different orbits can then give rise to the mass distribution 
we are interested in. 

Specifically, let A^orbits be the total number of orbits and 
A^coiis be the number of spatial cells (both will be specified 
below). For each orbit j £ A'orbits, we count the fraction 
of time, denoted by Oij, that it spends in each of the cells 
i G A'ceiia- The occupation number Wj for each orbit j is 
then determined by the following set of linear equations 



l,...,Afc, 



(4) 



E 

where Mi is the mass in cell i expected from the given mass 
distribution. We have checked the sum equals 
unity for each orbit. There are now various choices of how 
to actually solve Eq. [J] liner programming (Schwarzschild 
1979; 1982; 1993), Lucy's method (Lucy 1974; Statler 1987), 
maximum entropy methods (Richstone & Tremaine 1988; 
Statler 1991; Gebhardt et al. 2003), or least-square solvers 
(Lawson & Hanson 1974; Merritt & Fridman 1996; Capuzzo- 
Dolcetta et al. 2007). We chose the least-square method (cf. 
Wang et al. 2008). 

Because of the symmetry of the mass distribution spec- 
ified by Eq. (1), it is sufficient to only consider mass cells in 
the first octant in our analyses. Following Merri & Fridmann 
(1996), we divide the first octant into cells of equal masses, 
i.e. the first octant of each initial system is divided into 21 
sectors by 20 shells, where every sector has the same amount 
of mass inside. The first octant is further divided into three 
parts by the planes z — cx/a, y = bx/a and z — cy/b (see 
Fig. [1] upper left panel). Each of these three parts is sub- 
divided into 16 cells by the planes ay/bx — 1/5, 2/5, 2/3 
and az/cx = 1/5, 2/5, 2/3 (see Fig. [T] upper right panel). 
Therefore, the total number of cells in the first octant is 
16 X 3 X 21 — 1008. However, we only consider the innermost 
912 cells (the inner 19 sectors) when solving Eq.|4]since only 




1 r, r. 




Figure 1. The first octant is divided by planes z = cx/a, 
y = bx/a and z = cy/b (upper left panel) into three parts. 
Each part is subdivided by planes ay/bx = 1/5, 2/5, 2/3 and 
az/cx = 1/5, 2/5, 2/3 into 16 cells (upper right panel). In the 
lower panel, curve A is the circle of the minimal radius of 1:1 res- 
onant orbit at the energy Ek , curve B is the zero- velocity surface 
of the energy Bj,. There are 10 dotted lines x = ztanO divide the 
values of 6 from 2.25° to 87.75°. The 15 diamonds equally divide 
the radius into 16 parts. The diamonds are the initial positions 
from which the x — z orbits are launched. 



a few orbits supply densities in the outermost sector of grid 
cells; we simply discard the cells of the outermost sector. 

In a spherical system, the total energy and the three 
components of the angular momentum are integrals-of- 
motion. However, in a triaxial system only the energy re- 
mains constant (Merritt 1980; VaUuri & Merritt 1998; Pa- 
paphilippou & Laskar 1998). We therefore use the 7/8 order 
Runge-Kutta algorithm (Fehlberg 1968) for orbital calcula- 
tions, with 100 orbital times (see below Section |3TT| for the 
full integration time of each orbit. We followed Schwarzschild 
(1993) and Merritt & Fridman (1996) in assigning initial 
conditions from one of two sets of starting points (cf. also 
Wang et al. 2008): stationary orbits with zero initial veloc- 
ity, and orbits in the x — z plane with Vx = Vz = 0, and 
Vy — y^2{E — $) / in the first octant. Note that there 
is quite a large number of non-symmetric orbits that will 
lead to artifacts in the procedure if only being considered 
in the first octant. To circumvent this problem and to keep 
the symmetry of the system, we refiect the orbits from the 
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boundaries of each octant. Note that by our method the 
computational workload is not increased as would be the 
case when calculating eight octants. 

As mentioned above, there are 16 x 3 cells in each sec- 
tor. To obtain enough orbits for the library, we sub-divide 
the cells by the midplanes of the cells once again, i.e., the 
midlines of each grid as seen in the upper right panel of 
Fig. [1] these midlines equally divide the cell at x = az/c 
and X = ay/b. Thus we have 4 x 16 x 3 = 192 sub-cells in 
each sector. The central points on the outer shell surfaces of 
the sub-cells are the launching points. Hence there are 192 
stationary starting orbits in each sector. The total energy on 
the fcth sector is defined as the outer boundary shell E^, and 
for the stationary orbits inside the fcth sector this amounts 
to Ek ~ ^{x, y, z). For the orbits launched from the x — z 
plane the initial energy is E^ = '^{x, 0, z), as shown in Fig.[T] 
(lower panel). This figure further shows that the radius of 
the inner shell (marked as curve A) is the minimal radius 
of 1:1 resonant orbits (x;y), and the outer shell (marked as 
curve B) is the zero velocity surface. We define 10 lines sat- 
isfying x = zt&nO where 6 lies within the range 2.25° to 
87.75°. Along the radial direction, we equally divide the ra- 
dius between two boundaries into 16 parts with 15 points, 
where those 15 points are the initial positions for the orbits 
launched from the x ~ z plane. Hence there are 150 x ~ z 
plane starting orbits. 

In summary, we have 192 stationary and 150 x — z plane 
orbits in each sector amounting to a total of A^'orbits = 6840 
and we use Acoiis = 16 x 3 x 19 = 912 cells for the generation 
of our orbit library. The energy in each sector is a constant 
which equals the potential energy on the outer shell surface. 
As a result, the energies for the systems can be considered 
'quantized' with each system having 19 'energy levels'. In 
Fig. [2] there are some examples of asymmetric orbits: the 
upper four panels are stationary starting orbits, and the 
lower four panels are orbits launched from the x — z plane. 
In order to generate equilibrium Initial Conditions for the 
A-body simulations to be presented in Section O we need 
to symmetrize the orbits by using 'mirror particles'. 

Finally, in Fig. [3] we show the integration of mass as 
a function of energy for the model with a mass of 10^ Mq. 
There we find that more than half of the mass stems from 
loop orbits with chaotic orbits contributing more than 1/3 
of the mass; box orbits therefore do not play an important 
role in the model. We further like to note that the self- 
consistency of the model in MOND has been examined in 
Wang et al. (2008), and Antonov's third law was applied 
to check the stability of the models initially. However, it is 
unknown whether or not Antonov's third law is also valid in 
MOND so far. The most direct way to check for the stabil- 
ity and investigate the evolution N-body of the system is by 
means of N-body simulation to be elaborated upon in the 
following sub-section. 

2.4 ICs for N-body systems 

In order to study the stability and evolution of the systems, 
we need to convert the orbits into an N-body model. Accord- 
ing to Zhao (1996), the number of particles Uj on the jth 
orbit is proportional to the weight of the orbit, i.e., for an 
N particle system there are WjN particles on the jth orbit. 
Here we sample the particles on the jth orbit isochronously 





■0 



Figure 2. The asymmetric non-zero weights orbits in the orbit 
library. The upper four panels: stationary starting orbits and 
the lower four panels: x-z plane launched orbits. All the left 
panels are orbits projected on x-y plane and right panels are 
projected on x-z plane. 



at ti 



X (j + 0.5),i = 0,1,2, 



total integration time of the jth orbit. To this extend we in- 
terpolate the positions and velocities from the 6-dimensional 
output data of the Schwarzschild orbits. We generate 



(5) 



particles on the }th orbit and symmetrize the particles in 
phase-space. The remaining particles are kept as our Initial 
Conditions of the N-body system. 



3 N-BODY SIMULATIONS IN MOND 

All results presented in this section were obtained by evolv- 
ing our systems forward in time with using the N-body 
particle-mesh code NMODY (Ciotti et al. 2006; Nipoti et 
al. 2007a). 



3.1 Technical Details 

In our simulations, we have A = 8 x 10^ particles for 
each model, and choose a grid for the numerical integra- 
tion of Eq. [2] with 64 x 32 x 64 cells in the spherical co- 
ordinates {r,6,%l)), where the radial grids are defined by 
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Table 1. Total simulation times = 200 X Tsin 
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Figure 3. The accumulated of energy distribution of different 
orbit families in the intermediate model with a total mass of 
M = 10" M0. The horizontal axis is dimensionless energy of 
GM/i'ofcpc ' ^^"^ vertical axis is the integration of mass as a func- 
tion of energy. The dashed, dotted, dot-dashed lines are for box, 
chaotic and loop orbits. The solid line is for all of the orbits. 



Model 


N-body run duration T 


unit time Tgimu 


lO^O Mq 


0.94 Gyrs 


4.7 Myrs 


lO'^M© 


3.0 Gyrs 


14.9 Myrs 


10* Mq 


9.4 Gyrs 


47.0 Myrs 




Ti = 2.0tan[(i -f- 0.5)0.57r/(256 + l)]kpc. The density is ob- 
tained via a quadratic particle-mesh interpolation and the 
time integration is performed by the classical two order leap- 
frog scheme. As our time unit for all subsequent plots we use 
the following definition (cf. Wang et al. 2008) 

=( — ) =^-'^''y^[wnr,) [ikfc) 

which represents the Newtonian (or Keplerian) dynamical 
time at the radius of r = a without the factor of 27r. We 
remind the reader that the parameter Tsimu is neither the 
dynamical time nor the orbital time in general MOND sim- 
ulations. The orbital time in our MONDian systems is de- 
fined as the period of the 1:1 resonant orbit in the x — y plane 
(Wang et al. 2008). Fig.Ushows the periods of circular orbits 
at different radii for the three models presented here. Fig. |3] 
as well as Eq. |6] imply that the MONDian dynamical time 
at the radius of 1 kpc is about 7Tsimu, STsimu and S.STsimu 
for the three models whose masses are 10^° Mq, 10" Mq and 
10* Mq. We further like to note that the internal time step 
used by the code NMODY to integrate the equations-of- 
motion is , , where the factor 0.3 is a typical num- 

.^max |V-g| 

ber used in N-body simulations, and V-g = A-nGpeS, where 
Pe// is the effective dynamical density of the system, i.e. the 
sum of the baryon and (phantom) dark matter density in 
the Newtonian force law to produce the gravity or potential 
of baryons in MOND (see We et al. 2008). The time steps 
here are determined by the maximum values of V • g, which 
means the densest dynamical region of the models, where 
gravity changes most sharply. Note that all particles share 
a common time step that typically is 0.005 ~ O.OSTsimu. 

A flowchart of the technical steps involved in the process 
prior to the analysis stage can be viewed in Fig. [5] This figure 
summarizes the methodology of how to generate and evolve 



Figure 4. The period of 1:1 resonant circular orbits on the x — 
y plane as a function of radius. The solid, dotted, and dashed 
lines are for models with mass of IO^^Mq, 10"Mq and 10* M©, 
respectively. 

(6|)he A''-body systems. In Table [1] we present the total times 
each systems has been evolved for. 

3.2 Virial Theorem 

The scalar Virial theorem, W + 2K = 0, is valid for systems 
in equilibrium, where W is the Clausius integral, 

Vf/ = y" px-V^d^x, (7) 

and K is the kinetic energy of the system (Binney & 
Tremaine 1987). In the left panel of Fig.|S] we show that the 
evolution of —2K/W for all models is always about unity, 
as expected for an equilibrium system. We though note that 
during the first circa five Keplerian dynamical times (1.0 
kpc) all systems are moving from a quasi-equilibrium state 
with ~2K/W ^ 1.1-1.2 to -2K/W = l.OiO.l afterwards 
(marginally oscillating about unity). This figure demon- 
strates that our A''-body ICs start off in quasi-equilibrium 
and after approximately five Keplerian dynamical times (1.0 
kpc) can be considered fully relaxed. These 'hot' A^-body ICs 
could be due to a number of reasons including the resolution 
of the simulation and chaotic orbits, respectively. Regarding 
the latter, we need to mention that we compute the orbit 
library for 100 orbital times, and the time integration may 
not be long enough to ensure a relaxation of those chaotic 
orbits; particles coming from the chaotic orbits could lead 
to higher pressure 'overheating' the system. 
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Density (Eq. I) 



E 



(8) 



MOND Poisson 
Equation (Eq. 2) 
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Figure 5. Flowchart of the simulations. 




Figure 6. The evolution of 2K/\W\ for all three systems. The 
solid, dotted, and dashed lines are for models with mass of 
1O^°M0, IO^Mq and IO^A/q, respectively. The evolution is 
shown for 200 Keplerian dynamical times (1.0 kpc). 



In Fig. |6l we plot the velocity dispersion Wrms for all 
systems as a function of the simulation time unit. The plot 
indicates that each Wrms decreases by about 10% during the 
relaxation process and stays constant afterwards (with tiny 
variations though) Q 

Note that (as inferred from the right panel of Fig. [6]) 
the kinetic energy of the systems decrease nearly one quarter 
for the maximal evolved case after the relaxation0 However, 
the Virial theorem is still satisfied. That does not mean the 
energy conservation law is broken: W + K is not the total 
energy of a MONDian system, and it isnot conserved ei- 
ther. The total energy is still the conserved quantity but for 
a MONDian system it is given by (Bekenstein & Milgrom 
1984): 



^ There are typos in Wang et al. (2008) about the total mass of 
models and 

^ The kinetic energy K is proportional to v'^^^. 



where L is the Lagrangian of the MONDian system, defined 
by 



(V$)' 



(9) 



and ^{x^) is an arbitrary function with ii{x) — !F'{x'^). For 
an isolated system in MOND, the potential is logarithmic 
thus the potential energy is infinite. Therefore, the only 
meaningful quantity is the difference in energies between 
different systems (Bekenstein & Milgrom, 1984; Nipoti et 
al. 2007). However, the evident evolution oi W + K aX the 
very beginning (i.e. the first 5 Keplerian dynamical times 
(1.0 kpc)) shows that the N-body ICs are not accurately in 
equilibrium, and hence referred to as quasi-equilibrium. 

3.3 Energy Distribution 

One of the characteristic quantities to describe relaxation 
processes is the so-called differential energy distribution, 
i.e. the quotient of mass dM over the energy band inter- 
val \E,E + dE] (Binney & Tremaine 1987). The energy of 
a unit mass element is E = + (f>{x), where (j) is log- 
arithmically infinite in MOND and hence all particles are 
bound. But since the absolute value of potential energies is 
meaningless, we can define the zero point as the last point 
of the radial grid. Hence, there are positive relative energies 
for part of the particles though all of them are bound to the 
system. 

The left panels of Fig. [7] show the evolution of over 
200 Keplerian dynamical times (1.0 kpc) for all three mod- 
els (upper to lower) and we find that all distributions are 
rather similar. And the most pronounced evolution of the 
energy distribution is at the low-i5 end, where particles are 
most strongly bound to the system. All of the differential 
energy distributions have 19 peaks, as can be seen in left 
panels of Fig. [7] That is due to the energy definition of the 
Schwarzschild technique outlined in ij2.3l Inside every sector, 
the energy (kinetic plus potential energy) is a constant, while 
the outer shell is the zero- velocity surface of this sector. The 
adjacent two sectors have energy jumps at the shell. There- 
fore there are 19 'quantized energy levels' for our models, 
and for each model there are no mass distributions outside 
these 19 constant 'energy levels' and hence they appear as 
'valleys' in the left panels of Fig. [T] That explains why the 
curves appear noisy. We note that after the initial relaxation 
of about five Keplerian dynamical times (1.0 kpcjfjthe low- 
E end of the distribution becomes devoid of particles, i.e., 
particles are leaving the central regions where the potential 
well is deepest. This actually hints at a possible flattening 
of the initially present density cusp 7 = 1! We return to this 
issue later in sub- section 13.51 

Comparing the three left panels in Fig. [T] we observe 
that the system in the mild MOND regime (i.e. the model 

^ Even though we do not show the curves for 5 Keplerian dy- 
namical times (1.0 kpc) we acknowledge that the drop happens 
during that initial relaxation phase. Here we care about the long- 
term evolution within 200 Keplerian dynamical times (1.0 kpc) 
and hence decided to rather focus on the late evolution of the 
systems. 
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Figure 7. Left panels: Evolution of the diflferential energy dis- 
tribution . The panels (from upper to lower) correspond to our 
models with total mass of IO^^Mq, 10^ M© and 10* Mq, respec- 
tively. Right panels: The accumulation of energy distribution. 
The black lines denote the ICs, and the violet, blue, yellow and 
green lines show the differential energy after 50, 100, 150, 200 
Keplerian dynamical times (1.0 kpc). Both and the Energy 
are given in units where G=l and M=l. 



with a mass of 1 x 10^" A/©: upper-left panel) has the most 
significant evolution, whereas the model in deep MOND 
evolves least (lower- left panel). We therefore conclude that 
our ICs are most stable for the deep MOND regime. 

As seen in the right panels of Figure[7l the accumulation 
of energy distribution clearly confirms the previous conclu- 
sion from the differential distributions. After the relaxation, 
the mass in the inner region escapes to the outer, while the 
outer part is nearly unchanged. The mass distribution obvi- 
ously does not evolve after relaxation. 



3.4 Kinetics 

To further check upon the stability of our systems, we cal- 
culate the radial velocity dispersion profiles (Tr(r) as well as 
the anisotropy parameter 

*) = i-^^ (10) 

Here r is the spheroidal radius, the same as in Eq.[T] and ag, 
cr^ are the tangential and azimuthal velocity dispersions. 

The results can be viewed in the left panels of Fig.|8] We 
find that, for all three models, ar drops within the central 
2 kpc during the relaxation at the beginning of the simula- 
tion, i.e. the first five Keplerian dynamical times (1.0 kpc) 
(though not shown for clarity). Afterwards, there is only 
very little evolution noticeable. The reduction of ar in the 
core means that the ICs are too hot in the radial direc- 
tion to sustain equilibrium. We also note that the drop is 
more pronounced the less MONDian the ICs are. In the mild 
MONDian model with a mass of M = W^^Mq, (upper-most 
panel) the model obviously appears to be hotter inside than 
outside. This trend is weakened for the deep-MOND model 
with a total mass of M = 10* Mq where the self-gravity of 
the system is much weaker than oq. The slope of ar{r) os- 
cillates around a constant value of approximately 20 km/s. 
In the intermediate model with M = W^Mq the slope is 
between the most massive and least massive ones. 

However, after the relaxation, the slope of ar{r) has 
the same behavior in all three models, radially cooling down 
towards the cores. The curves of velocity dispersion after re- 
laxation look like the rotation curves of disc galaxies at the 
similar mass range (Milgrom & Sanders 2007; Gentile 2008). 
In the centres, the tangential velocity dispersion ag + 
plays an important role for the Virial Theorem and keeps 
the cores in equilibrium. Therefore, the systems prefer more 
isotropic velocity dispersions in the cuspy centres. To con- 
firm this, we also present the anisotropy parameter /3(r) in 
the right panels of Fig. [S] We do find the expected small 
values of f3 in the centres as well as a radial increase of 
/3(r). Inside 1 kpc the /3-profiles oscillate during the whole 
evolution while they remain stable in the outer parts. Fur- 
thermore, the anisotropy increases with radius in all three 
models out to about 25 kpc where it turns approximately 
constant, (3 = 0.6, i.e., the velocities are distributed hyper- 
radially. 

Nevertheless, within 2 kpc, there should be a more sub- 
stantial redistribution of kinetic energies in the tangential 
direction after relaxation due to the evolution of ar seen in 
the left panels of Fig. [S] otherwise, the systems lose quite a 
lot of kinetic energy in the core. Unfortunately, the evolu- 
tion of f3 in the central region, seen in right panels of Fig. [S] 
is not as large as expected. Thus, there are outfiows of ki- 
netic energy from the centres of the systems. The values of 
velocity dispersion and anisotropic parameters in the deep 
MOND regions of our systems fit pretty well with the ana- 
lytical predictions of isothermal spheres by Milgrom (1984; 
1994). 

3.5 Mass distributions 

Due to parts of the kinetic energies spilling out of the cores 
(cf. sub-section 13.31 and 13. 4p . the mass densities could redis- 
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Figure 8. Left panels: Evolution of the radial velocity disper- 
sion rjr{r). Right panels: Evolution of the velocity dispersion 
anisotropy /3{r). The upper, middle and lower panels are cor- 
responding to models of M = 1.0 X 10^°Mq, M = 1.0 X IO^Mq 
and h[ = 1.0 X IO^Mq. The ordering of the panels corresponds 
to Fig. [7] as does the colouring of the lines. 



tribute at the same time. Indeed, we find that there are out- 
flows of mass: the cuspy centers with an initial value of 7 = 1 
are flattened. This can be viewed in Fig.|9]where we show the 
densities along the major axis (left panels) and cumulative 
(right panels) mass distributions for our MOND models. 
With regards to the density panels, the three models show 
a similar behavior. It is clear that the mass is redistributed 
during the relaxation with losses in the very central region 
of r < 0.5 kpc and gains outside. The density curves are 
ocsillating around the initial analytical density as given by 
Eq. [1] Therefore, the system becomes slightly less cuspy and 
keeps the triaxial density after reaching equilibrium. Note 
that the density distribution still remains triaxial after the 
system is in equilibrium, and there is no obvious evolution 
within 200 Keplerian Dynamical times (at the typical scale 
a=1.0 kpc). The right panels of Figure[9]show the total mass 
inside the radial direction r. The black dashed straight lines 



in right panels are defined by Mo = the mass to pro- 

duce the gravity acceleration ao in a point mass approxima- 
tion. A'lo is the watershed of the enclosed mass producing 
MOND and Newton dominating gravities. At a certain ra- 
dius r, when the enclosed mass is smaller than Mo, there 
occurs a transition to MONDian gravity. We find that in 
all of the three models MONDian eflects cannot be ignored. 
Even for 1.0 x 10^° Mq, the MONDian gravity dominates 
the regions of r > lO"'^ ~ 2 kpc. Obviously, the model 
M = 1.0 X 10** M© is in deep MOND region. The colours 
show the evolution of the systems. Mass in the inner part of 
the system is lost during the density re-distribution, while 
in the outer part, beyond 4 kpc, the total mass is not af- 
fected. We further note that after the redistribution (during 
the relaxation) the mass distribution has stabilized. 



3.6 Shape 

As confirmed in H3.5I the mass redistributes inside our sys- 
tems during relaxation. Hence, the question arises whether 
the shape (i.e., the initial triaxiality) remains stable or un- 
dergoes changes. To address this, w e show t he evol ution of 
the axis ratios of the eigenvalues \fT^^^^jL^ and ^ Izz/ Ixx 
of the moment of inertia tensor rriijXiXj, {rriij = M/N) in 
Fig. 1101 The three models give similar results, hence we high- 
light the model with a total mass of W^Mq in the upper 
panel. We find that the axial ratios (as a function of radius) 
merely evolve about 10% within 200 Keplerian dynamical 
times (1.0 kpc). The system is rounder in the center, but 
the whole system keeps the triaxial shape during the long 
stage of evolution. It is clear that the axis ratio between the 
minor and major axes is more stable than that of the in- 
termediate and major axes. Not only the ratios are almost 
constant, but also the absolute values of Ixx, lyy and Izz 
seem in dynamic equilibrium and stable, changing less than 
20% during the oscillation (cf. Fig llip . However, we note 
that at the time t — the ratios of lyy/Ixx and Izz/Ixx 
do not accurately equal b : a = 0.86 and c : a — 0.7 in most 
of the inner regions, which is caused by the numerical effects 
in generating the A'^-body ICs. However, at the edge of the 
galaxy (i.e., including more than 80% of the total mass), the 
axis ratios are close to the suggested ones. We conclude that 
the ICs generated by our application of the Schwarzschild 
technique roughly lead to a 5% error for the axis ratios. 

A study of the tensor kinetic energies Kxx, Kyy and 
Kzz, defined as Kxx ~ 0.5 < Vx ■ Vx >, shows a similar 
behaviour to the moment of inertia tensor analysis presented 
above. The ratios remain constant although the absolute 
values change by at most 20%. This can again be verified in 
Figure [TT] where we plot the inertial (left panel) and kinetic 
energy (right panel) components for the three models. 

As a final note, considering existing Schwarzschild plus 
N-body simulations in the literature, we find that the evo- 
lution seen in our MONDian cuspy elliptical models is com- 
parable to that seen in Fig. 5 of Poon & Merritt (2004, ApJ 
606, 774) for triaxial ellipticals in Newtonian gravity. Our 
simulations are much longer than 10 crossing times, which 
provides a typical scale for checking stability in Newtonian 
Schwarzschild simulations in the literature (Poon & Merritt 
2004, Zhao 1996). 



N-hody simulations for testing the stability of triaxial galaxies in MOND 9 




0,1 1,0 10,0 -1,0 -0,5 0,0 0,5 1,0 1,5 2,0 

< (kpc) log,„ rodius (kpc) 



Figure 9. The evolution of the mass distribution for the models 
with M = lOl°M0 (upper), M = 10^ Mq (middle), M = 10* Mq 
(lower). The left panels show the density distributions on the 
major axis the density information can be obtained from the axis 
ratios of Figure FlOl The right panels show the accumulated mass 
inside the radius r. The dashed black lines in the right panels 

2 

are defined as , which are the watersheds of enclosed mass 
producing MONDian dominating gravity (below the lines) and 
Newtonian dominating gravity (upon the lines). The colouring of 
the lines is representative of the evolutionary stage of the model 
and corresponds to Fig. [7] 



4 CONCLUSIONS AND DISCUSSION 

We explored the stability and evolution of the triaxial 
Dehnen model (Dehnen 1993; Merritt & Fridman 1996; 
Capuzzo-Dolcetta et al. 2007) with a 7 = 1 central 
cusp using MOND. We utilized the Schwarzschild method 
(Schwarzschild 1979) to build orbit models which were in 
turn used to generate initial conditions (ICs) for N-body 
simulations using the method outlined in Zhao (1996). These 
ICs were evolved forward in time for 200 Keplerian dynam- 
ical times (at the typical length scale of 1.0 kpc) by the nu- 
merical integrator NMODY developed by the Bologna group 
(Ciotti et al. 2006; Nipoti et al. 2007) and designed to in- 
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Figure 10. Evolution of axis ratios with the median model of a 
total mass of 1.0 X 10^ M0 (Upper panel), 1.0 X W^°Mq (lower 
left panel) and 1.0 X 10* (lower right panel). The lower and 
upper series of lines are for the ratios of minor : major axis and 
intermediate : major axis, i.e., Izz/Ixx and ^ lyy/Ixx- The 
different line symbols are defined the same as in the figure: solid, 
dotted, short dashed, dot-dashed and long dashed lines are for 
system evolving 0, 50, 100, 150 and 200 imu^ where Tg^jy^u is 
the Newtonian orbital time at 1.0 kpc . 



elude the effects of MOND. We additionally ran the same 
simulations with a second MONDian gravity solver AMIGA 
(Llinares, Knebe & Zhao 2008, cf. Appendix[B]| based upon 
an entirely different grid-geometry to confirm the credibility 
of our results. 

In our simulations, the virial theorem was satisfied 
at all times. We showed that the systems start in quasi- 
equilibrium with a short relaxation phase of approximately 
less than five Keplerian dynamical times (1.0 kpc). We found 
outflows of energy and mass from the centres of the sys- 
tems under investigation. Hence, during the relaxation stage, 
there is a flattening of the initially present 7=1 cusp to 
a core. Despite the obvious mass redistribution, we need to 
acknowledge that the shape of the systems remained un- 
changed in the course of the simulations; the axis ratios of 
the eigenvalues of the moment of inertia tensor (as well as 
the kinetic energy tensor) stayed constant. 



10 Wu et al. 




.. i .... i .... i ... . lOOOL 
50 100 150 200 50 100 150 200 

Keplerian dynamical times at 1 .Okpc Keplerian dynamical times at 1 .Okpc 




50 100 150 200 50 100 150 200 

Keplerian dynomicol times at 1 .Okpc Keplerion dynamical times ot 1 .Okpc 





80 






200 






















60 


r~\ /\ l„ /~\ _ 




150 


















(kpc 






1 








o 
















40 






100 


















D 






























20 






50 





































50 100 150 200 

Keplerian dynamical times at 1 .Okpc 



50 100 150 200 

Keplerian dynamical times at 1 .Okpc 



Figure 11. Upper, middle and lower panels are different mass 
models of IO^^Mq, lif Mq and 10** M©. The left panels are 
the evolution of the inertial tensor Ixx (solid line), lyy (dotted) 
and Izz (dashed). The right panels are the evolution of kinetics 
energy K^x (solid), Kyy (dotted) and K^z (dashed). The total 
simulation time is 200 Tgimu ■ 



The effects of resolution of the simulations should not 
remain unmentioned. We found that the potential calculated 
from the N-body ICs differs by 10% compared to the analyt- 
ical potential. Furthermore, the analytically predicted veloc- 
ity dispersions at the initial time are 107.3 km/s, 54.2 km/s 
and 29.3 km/s for the models M = 10^", lO'', 10* Mq, respec- 
tively. However, they do not match the (numerical) values 
plotted in right panel of Fig. [B] Hence we use the Clausius 
integral \ W\ in Equation [T] calculate the analytical densi- 
ties for the systems at t = 0, to minimize the errors. More- 
over, we have found that due to the resolution limitation of 
the NMODY code, the errors accumulate during the sim- 
ulations, which is insensitive to more massive systems, but 
becomes an issue when the mass of the system decreases. 
This causes small non-zero net velocities in the systems. The 
simulation centres of the systems with W^Mq and 10* 
move significantly after 200 Keplerian dynamical times (1.0 
kpc) and hence we restrict our analysis to this time frame. 



To further check the credibility of our results and the de- 
pendence on the code, we ran the simulations again with 
a technically substantially different code (AMIGA), which is 
also capable of integrating the analog to Poisson's equation 
(cf. equation [2] in Appendix iB]! . The results are practically 
indistinguishable reassuring their tenability. 

We like to close with the notation that our systems are 
isolated systems, corresponding to the cases of field galaxies. 
The self-potentials of the systems in MOND are logarithmic 
at large radii, therefore no stars can escape from such sys- 
tems. However, for any system embedded in external fields, 
the potential is truncated when the strength of the exter- 
nal field becomes comparable to the internal field (Milgrom 
1984; Wu et al. 2007). Therefore, Poisson's equation should 
be modified to 



|V$.. 



- gext\ 



ao 



(V$i„t - gcxt) 



= 4TTGph, (11) 



where the /.t-function is determined by both the internal 
and external gravitational accelerations. Hence the strong 
equivalence principle is violated, and the directions along 
and against the external field, the /i-function has different 
values even though the mass density distributions are the 
same. A direct result is that the potentials become non- 
symmetric along and against the directions of the external 
field, i.e., a symmetric system is not in equilibrium due to the 
non-symmetry of self-potential. Therefore, MOND predicts 
that there are no real symmetric systems within the external 
gravity backgrounds. This will be explored in greater detail 
in a future paper (Wang et al. in preparation). 
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APPENDIX A: SYMMETRY AND 
NUMERICAL CHALLENGES 

Evolving the systems without filtering the high frequent 
components of mass during the Legendre transformation in 
time for up to 200 Keplerian dynamical times (1.0 kpc) we 
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Figure Al. The average values of positions (z) vs. velocities (v^) 
projected on the z axis for each non-zero weight orbit. The model 
has the total mass of 10^ Mq. 



Figure A2. Left panel: The shifted positions of the centre of 
mass of the three systems. The solid, dotted, dashed lines are 
for the systems with their total mass of IO^^Mq, 10® Mq and 
10** Mq. The colours of black, yellow and blue correspond to the 
z, y and x axis. The right panel: net velocities of the systems in 
the three axis directions. 



find that they appear 'unstable'. However, a detailed inves- 
tigation revealed that is a numerical effect rather than a 
physical instability. There is a small, uneven force in the 
z direction which comes from the asymmetry of the den- 
sity distribution of ICs for N-body and systems during the 
simulations and the errors obviously accumulate when the 
numbers of particles are not big enough. Further, the total 
momentum of the systems is not conserved giving a non-zero 
net velocity along the minor axis. Hence the code requires 
a large number of particles for smooth density distributions 
to make sure the tiny asymmetry of some particles does not 
affect the whole system. 

Note that this effect is more serious when the systems 
are not symmetrized by utilizing 'mirror particles' inside the 
systems. It is known that particles generated from asymmet- 
ric orbits (e.g. the 'banana orbits') could break the symme- 
try of the systems in phase space. Furthermore, there are 
a couple of hundred of chaotic orbits with positive weights, 
and, during the Schwarzschild process, the time integration 
of 100 orbital times may not be long enough to obtain sym- 
metry. We therefore show in Fig. lAll the average value of po- 
sitions and velocities along the z-axis for every orbit of the 
A'^-body ICs of the model with mass l(f Mq. We plot z vs. Vz 
since the z-direction displays the most serious shifting. For a 
perfectly symmetric system the values of z and should be 
close to zero, while we find that they are not. Hence we need 
to symmetrize the systems prior to the N-body procedure. 
The simplest way to achieve this is by placing 'mirror parti- 
cles' into the system, i.e. using a minus sign in front of the 
6-dimensional components. Therefore, the total numbers of 
particles increases to iV X 2^* = 64iV. 

To illustrate (and quantify) this effect we plot in Fig. lA2l 
the centre-shifts along the three axes (left panel) as well as 
the evolution of the net velocities along the same axes (right 
panel) during the first 40 Keplerian dynamical times (1.0 
kpc) for the models without the 64 mirrors. We find that the 
most massive system (i.e. lO^^M©, which is in mild-MOND 
gravity) is least affected by these numerical artifacts. As a 
matter of fact, this particular system shows credible signs 
of stability even after 200 Keplerian dynamical times (1.0 
kpc). We need to acknowledge that this is partly due to 
our definition of the time unit (cf. Eq. [6)l : it is shorter for 



more massive systems. However, all the analysis presented 
in this paper indicates that the system is stable despite the 
apparent numerical artifacts of the code. We though cannot 
evolve the system further in time as the accumulation of 
errors would lead to substantial deviations from the system's 
equilibrium state; but this 'instability' is caused by numerics 
rather than physics! Given the technical particulars of the 
NMODY code such a centre-shift will lead to a decrease in 
resolution as the code utilizes a spherical grid. 



APPENDIX B: COMPARISON WITH 
ANOTHER MOND SOLVER 

As just highlighted in Section [X] there are numerical chal- 
lenges to evolving our systems under MONDian gravity us- 
ing the N-body code NMODY. In order to confirm that the 
results are not unique to this one code we therefore decided 
to also use another novel solver for the MONDian analog to 
Poisson's equation, namely the AMIGA code (Llinares, Knebe 
& Zhao 2008). AMIGA is the successor to MLAPM (Knebe, 
Green & Binney 2001) that has recently been adapted to also 
solve Eq. I^PI The code utilizes adaptive meshes in Cartesian 
coordinates in a cubical volume as opposed to the spherical 
grid of NMODY. The solution is obtained via multi-grid re- 
laxation and we refer the interested reader to Llinares et al. 
(2008) for more details. 

However, here we need to elaborate upon the boundary 
constraints as we cannot assume that the potential on the 
boundary will be a constant: the box is a cube and not a 
sphere. We decided to use the solution for a point mass in 
the center of the box 

*W = -^^(---) +(/W-/M), (Bi) 

2 \r roJ 

with 



* We like to note in passing that MLAPM has already been suc- 
cessfully applied to study cosmological structure formation under 
MOND (Knebe & Gibson 2004) under certain assumptions. 
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fir) 



ao 



q = 



— +4r2 +ln 
2r 



GM 
ao ' 



4r' 



(B2) 



where, M is the total mass in the box and ro is a length 
scale (a constant of integration). For ao we recover 
the Newtonian solution and for ao finite and r — > oo we 
have ln(r), which is the typical behaviour for any MONDian 
solution. In the case that we use ro = B, with B being the 
size of the cubical box, we end up with "1> = in a sphere 
of radius B, that is equivalent to the conditions used in 
NMODY. 

We now run simulations with the same Initial Con- 
ditions for N-body as used with NMODY utilizing a do- 
main grid with 128'^ cells. Each of these domain grid cells 
is refined and split into eight sub-cells once the number of 
particles inside that cell is in excess of 6. The box size is 
B = 165.5152kpc and the scale for the boundary conditions 
is ro = 82.7576 Kpc (half of the box). 

The results obtained are similar to the NMODY simula- 
tions. The system is stable with a normal secular evolution. 
We observe the same kind of evolution. We confirm that all 
other quantities behave in a similar manner too, and hence 
are confident that the results presented in the previous sec- 
tion[3]are not dominated and/or contaminated by numerical 
artifacts. 



APPENDIX C: LONGER TIME EVOLUTION 
OF THE 10^°Mq model 

We initially ran the models for 200 simulation time units 
(i.e. 200 circular orbital times at the length of 1.0 kpc, see 
Table [1} • Hence the most massive model has been evolved 
for the least time, about 1 Gyr. For brighter galaxies because 
they are also bigger galaxies, the dynamical time is longer. 
This is the opposite to our trend in Table 1. This can be 
understood because our toy galaxies do not sit on the funda- 
mental plane. To see this, we estimate r^, the characteristic 
length of an elliptical galaxy on the fundamental plane (e.g., 
Eq. 9 in Zhao, Xu, Dobbs 2008, which follows Faber et al. 
1997 ). 



log 



GMr- 



350 X 10-iOm/s^ 



= -1.52 log 



M 



1.5 X IO^Mq 



±0.5, (CI) 



we find r-h = 0.082+^;^*° kpc for the model with 10'" Mg. 
This is one order of magnitude smaller than our assumed 
size 1 kpc. The dynamical time for a 10^°Mq galaxy sitting 
on the fundamental plane is about 0.1 Million years. 

To be on the safe side, we re-ran the model for about 
3 Gyrs {650Tsimu), and show its virial ratio (Fig. IClj ). And 
our conclusion in the i]3.2l does not change. The virial ratio 
oscillates around 1 within 10% at most. 
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